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ABSTRACT 

Demanding aerodynamic modelling requirements for military and civilian aircraft have motivated 
researchers to improve computational and experimental techniques and to pursue closer collaboration in 
these areas. Model identification and validation techniques are key components for this research. This paper 
presents mathematical model structures and identification techniques that have been used successfully to 
model more general aerodynamic behaviours in single-degree-of-freedom dynamic testing. Model 
parameters , characterizing aerodynamic properties , are estimated using linear and nonlinear regression 
methods in both time and frequency domains. Steps in identification including model structure determination, 
parameter estimation, and model validation, are addressed in this paper with examples using data from one- 
degree-of -freedom dynamic wind tunnel and water tunnel experiments. These techniques offer a methodology 
for expanding the utility of computational methods in application to flight dynamics, stability, and control 
problems. Since flight test is not always an option for early model validation, time history comparisons are 
commonly made between computational and experimental results and model adequacy is inferred by 
corroborating results. An extension is offered to this conventional approach where more general model 
parameter estimates and their standard errors are compared. 


NOMENCLATURE 
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a , bi, c 
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Ci 

C N 


c c 


transfer function coefficients 
Fourier coefficients 
indicial function parameters 
wing span, ft 

rolling-moment coefficient 
normal-force coefficient 
mean aerodynamic chord, ft 
in-phase components of C a 


D 

do, dj 

f 


out-of-phase components of C a 

time domain differential operator 

constants 

frequency, Hz 


m 

k 


m 

N 


p> q 

R 2 


s 

t 

V 

y 

z 


sum of squared residuals 

reduced frequency, 2 nfi / V 

half chord or half span, ft 

No. of harmonics in Fourier expansion 

number of data points 

roll and pitch rates, rad/sec 

squared multiple correlation coefficient 

estimated standard error 

time, sec 

airspeed, ft/sec 

output variable 

measured variable 
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a 

= angle of attack, rad or deg 

T 

= dummy integration variable 

OC 0 

= mean angle of attack, rad or deg 


= non-dimensional time constant, — 

V 

p 

= sideslip angle, rad or deg 

9 

= state variable 

V 

= residuals 

9 

cP 

= vector of unknown parameters 
= variance 

CO 

= angular frequency, rad/sec, 2nf 


Subscripts 

A = amplitude 

a = N or / 


Superscripts 

= complex variable 
= complex conjugate 


Aerodynamic Derivatives 

S (oo) ^ Cl P = ~d^' Cn » (oo) s Cn <x = 
c h (co) ^ Ci p = V^’ CiV ? (oo) s 


da 


1.0 INTRODUCTION 

S ince the early days of flight, aerodynamicists have focused on the problem of finding adequate 
aerodynamic models with good prediction capability. The problem has become more acute over time with 
modern military aircraft requirements expanding the flight envelope, demanding more capability in nonlinear 
unsteady flight regimes, and civilian aircraft requirements for safety, demanding more capability in similar 
adverse aerodynamic conditions. Over time researchers have endeavored along two paths to address this 
nonlinear unsteady modelling problem: a numerical path using high fidelity Computational Fluid Dynamics 
(CFD) technology and an experimental path using direct measurements of aircraft responses in flight or scale 
models in wind tunnels along with System Identification (SID) technology to extract adequate mathematical 
models from the measured data. Advances in both numerical and experimental technologies have created 
opportunities for the combination of these technologies to make significant improvements in handling more 
difficult aircraft modelling problems. 

In order to capitalize on these opportunities, model validation and identification should be carefully 
considered. Although it is not common practice for wind tunnel or CFD results to be presented with associated 
error bounds for validation comparisons, some efforts, Refs. [1-4], are being made to improve this 
shortcoming. Model identification, considered in this paper, is a research area that involves a number of 
disciplines in modelling and measurement science. The subject of aircraft system identification, in particular, 
is documented in Ref. [5] and MATLAB® software is included for a variety of system identification methods 
commonly used in aerospace engineering practice. The identification process includes a number of key steps 
such as, experiment design, model structure determination, parameter estimation, and model validation. These 
steps lead to specific demands on test techniques and test facility capabilities. The demand for more general, 
higher fidelity, models has highlighted some limitations of the conventional approaches. For example, in 
ground-based dynamic wind-tunnel experiments where single frequency sinusoidal inputs are typically used, 
wide-band inputs have been shown to be more effective and less costly since fewer runs are required, Ref [6- 
7]. Advancement in these areas is hampered by the expense of developing advanced test facilities capable of 
more general motions. CFD offers an opportunity to ameliorate a number of these limitations. 

As an initial step toward improved integration of CFD and SID technologies, this paper presents identification 
methods used at NASA Langley Research Centre that can be applied when nonlinear unsteady aerodynamic 
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responses are present. The model identification problem for conventional one-degree-of-freedom motion in 
wind tunnel forced-oscillation experiments is described. Mathematical model structures are presented that 
include conventional static and rotary dynamic terms but are extended by including indicial functions to 
represent unsteady responses. Although other model structures have been studied, the authors’ goal is to 
maintain the conventional stability and control derivative structure to capitalize on the large engineering 
knowledge base built around that structure. Model parameters describing aerodynamic properties of the 
system under test can be estimated using a least-squares principle. The resulting estimators form linear or 
nonlinear regression problems that can be formulated in the time or frequency domain. Five parameter 
estimation methods are described and demonstrated on experimental data from wind tunnel and water tunnel 
experiments. These methods are harmonic analysis, nonlinear regression, two-step linear regression, equation 
error, and output error. Model validation of these techniques is demonstrated with examples using the NASA 
Generic Transport Model (GTM), Fig. la, and two scale models of the F-16XL: a 2.5% scale water tunnel 
model, Fig. lb, and an 18% scale wind tunnel model not shown. 


2.0 MATHEMATICAL MODEL STRUCTURES 

Estimation of aerodynamic parameters from wind tunnel data requires that a mathematical model of the 
aircraft is postulated. The mathematical model includes both the aircraft equations of motion and the 
equations for aerodynamic forces and moments, known as the aerodynamic model equations. For most of the 
nominal operating envelope, where the aerodynamic flows are attached and quasi-steady in behaviour, 
aerodynamic model equations are usually developed using time -invariant linear terms (stability and control 
derivatives). In off-nominal cases these terms are generalized by including time -dependent and nonlinear 
terms. Results from numerous studies, such as Ref. [6], have shown the dependence of aerodynamic 
parameters on frequency. This dependency contradicts the basic assumption that stability derivatives are 
constants. Because the effect of frequency is related to unsteady aerodynamics, the aerodynamic model 
equations will be formulated in terms of indicial functions or “step-response functions”, Refs. [8-9]. As an 
example, the normal force coefficient is considered as 


t £ t 

C N (0 = C N (0) + J C Na (t - z)d{z)dz + — J C Nq (t - z)q(z)dz 
0 2V 0 

( 1 ) 

where C N (t) and C N (t) are the indicial functions. For further development of Eq. (1) it will be assumed 
that 

a) only the increments to steady conditions are considered 

b) the effect of q(t) on the normal force can be neglected 

c) the indicial functions can take the form of a simple exponential 

C Na (t) = a(l-e- b ' r ) + c 

= C Na (a-j)-ae^ h]t 

where the indicial function is separated into a steady and deficiency function, Ref. [10]. Consequently, this 
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model can be simplified as 

C N (t) = C Nr/ (oo)a(t) + ^-C Nn (°0)q(t)-a\e-W- T) d(T)dT 
2V V Q 


or in operator form as 


(2a) 


£ CL 

C N (t) = C Na (oo )a(t) +—C N (°o)q(t) - Da 

V ^ J-J + D 'i 


By introducing 


(2b) 


r/(t) = \e ^ T ^a{z)dz 
o 


the state-space form of Eq. (2) is 


( 3 ) 


7j(t) = -b ] ?j(t) + d(t) 

£ 

C N it) = C N (oo )a(t) + —C N (oo )q(t) - ar/(t ) 

a V q (4) 

Applying the Laplace transform to Eq. (4), the transfer function for the coefficient is obtained as 


C yy is) As 2 + Bs + C 

ais) s + b x 


where s is the Laplace transform parameter and 


( 5 ) 


A = -C n (oo) 

v <? 

B = C N a (°°) - a + b \ ^ C N q (°°) 

C = b x C Na { co) (6) 

f is the half chord or half span as required. A different model structure suitable for the analysis of oscillatory 
data, with amplitude a A , follows from the steady solution of Eq. (2) as 


C N it) = a A C Na sin icot) + a A kC N ^ cos icot) 

( 7 ) 
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where k = 7rcf /V is the reduced frequency. For the conventional model form with linear aerodynamics, 
C N (t) is assumed to be functions of a,a,q , and q . The in-phase and out-of-phase components, C Na and 

C N , can then be expressed in terms of steady flow and acceleration terms 


Cyu — C N (oo) — k C 


Nn 


(8a) 


C N q -C Nq (co) + C N . 


(8b) 


These acceleration terms are the conventional equivalent to the unsteady terms. In order to explain the 
variation of the in-phase and out-of-phase components with frequency these components were expressed in 
terms of deficiency functions, as in Eq. (2). In Refs. [10-1 1], the model for these components has the form 


Qwv — (°°) — a 


z 2 k 2 


Na '~' "l + z 2 k 2 


(9a) 


C 


N n 


c 


N q (°°) ' 


l + r 2 k 2 


(9b) 


where z x = 2V / cb { is an aerodynamic time constant. In the aerodynamic model equations presented here there 
are six unknown parameters: C Na , C N , a , b x (or r^, C Na , and C N . Techniques for estimating these 
unknown parameters are provided in the next section. 


3.0 ESTIMATION METHODS 

The parameters in the aerodynamic model equations can be estimated in various ways, e.g., U.S. Air Force 
DATCOM, based on an extensive set of rules for aerodynamic dependencies and aircraft geometry, can 
produce basic aerodynamic parameters but these are limited to basic geometries and steady-flow dynamics. 
Alternatively, analytic methods such as CFD, where general physics-based equations define the mathematical 
model, can provide general aerodynamic solutions but application to stability and control problems is still an 
area of research. Although these methods work well to varying degrees, in particular for low angle of attack 
and low rotational rate cases, the best aerodynamic predictions are still obtained using experimental methods. 
Flight experiments provide the most direct measurements of aerodynamic behaviours but this is the most 
expensive approach and not available in early design phases. The wind tunnel provides the next best approach 
with current technology; however, this approach has a number of difficulties related to low speeds, similitude 
scaling, and tunnel effects. In this paper the aerodynamic parameters will be estimated from wind tunnel 
dynamic tests to demonstrate modeling methodology that can also be applied to appropriate CFD simulations. 
Five methods for estimation of parameters in aerodynamic model equations are presented in this paper: 
Harmonic Analysis, Nonlinear Regression, Two-Step Finear Regression, Equation Error, and Output Error. 
They are based on application of a least-squares principle to experimental data assuming that the aerodynamic 
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model structure is known. Wind tunnel dynamic tests usually include forced oscillations about the aircraft 
body axes at selected angles of attack, sideslip, frequencies, amplitudes, and Reynolds numbers. During these 
tests the inputs (Euler angles) and the outputs (aerodynamic coefficients of forces and moments) are recorded. 
All the estimation methods presented are demonstrated using data from wind tunnel experiments using the 
NASA GTM and 18% scale F-16XL, and from water tunnel experiments using a 2.5% scale F-16XL model. 

3.1 Harmonic Analysis 

A method of harmonic analysis, Ref. [12], is applied to measured aerodynamic coefficients to allow 
estimation of the in-phase and out-of-phase components. For the development of this method it is assumed 
that a periodic function y(0 = y(t + 2 n) with the period 2k is represented by a series of discrete values y(i) at 

, .. 2m , . _ 

t(i) = ,* = 1,2,..., A 

N 

It is assumed that y(i) can be approximated by a trigonometric series 


m m 

y{ 0 = 4 , j cos jco {) i + ^ B j sin jeo Q i 

j = i j = i 

( 10 ) 

2k 

where a> 0 = — . This series uses an orthogonal-basis of harmonic sinusoids. It is further assumed that the 
measurements of y(i) are obtained as 


z(i) = y(i) + v(i) 


( 11 ) 

where z(i ) are the measured values and v(i) is white measurement noise with zero mean and constant variance 
cr 2 . The parameters Aq , Aj and Bj , are Fourier coefficients that can be estimated from the measurements by 
minimizing the least squares (FS) criterion 


^L 5 (^)=S[z(0-j(0] 2 

i=l 


( 12 ) 

where 0 is the vector of unknown parameters, (A„, A h B h A m , B m ). Then the LS estimates of parameters in 
(10) are 
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2 N 

Aj = — J^z(i)cosjco 0 i 
2 N 

= — ^z{i)smja> Q i 


The parameter variances of these estimates are 


N 

: (*,H 2 (^H 2 2 


N 


for all j. The estimate of variance <x 2 is 


( 13 ) 


( 14 ) 


s 2 =Tz[z(0-y(0] 


where y(i) follows from Eq. (10) by replacing the parameters by their estimates. 


( 15 ) 


The adequacy of the model given by Eq. (10) can be assessed by the coefficient of determination that indicates 
how much variation in the data is explained by the model. R 2 is expressed as 


£[z(0-y(0] 2 

R 2 = 1- — , o < f? 2 < 1 

z[*(0 -z? 

i = 1 

where z is the mean value of z(/) . 


( 16 ) 


For the pitch oscillation case, considering only the first harmonic and a model with linear aerodynamics, the 
in-phase and out-of-phase components of C N (t ), can be expressed in terms of coefficients Aj and B u 
respectively, as 


^ 

C N = — and C 


a a 


N n 


A_ 

ka A 


( 17 ) 

where k = 7rcf /V is the reduced frequency. Eq. (17) can be easily modified for the roll and yaw oscillation 
case, Ref. [12]. 

An example of harmonic analysis performed on pitch oscillatory data is given in Fig. 2 for the NASA Generic 
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Transport Model (GTM). Ordinate values were removed in order to maintain proprietary agreements. In this 
figure the in-phase and out-of-phase components, and the coefficient of determination are plotted against 
angle of attack for seven frequencies and one selected amplitude. This experiment was conducted in the 
NASA Langley 14x22 Wind Tunnel. Model geometry for the GTM is given in Fig. la. 


3.2 Nonlinear Regression 

Harmonic analysis results shown in Fig. 2, indicate that the in-phase component, C Na , does not vary with 

frequency and, although not shown, its values agree with C Na determined from steady measurements. On the 

other hand, the out-of-phase component, C N ^ , shows variation with frequency for /< 0.43 Hz and 10° < a 0 < 

50°. Values of R 2 indicate that a linear aerodynamic model may be adequate for a 0 < 35°. In this case, three 
unknown parameters, C N ^ , a , and , in Eq. (9b), can be estimated by a nonlinear regression method, 

suggested in Ref. [13], using the frequency dependent measurements. 

A general model for this method can be formed as 


z(j) = g[x(j), 6] + v(j) j = 1 , 2 , m 


(18) 

where x(j) is a vector of regressors computed from measured data at the jth data point, g is a nonlinear 
function of x(j) , and unknown parameters are given by the vector, 6 . The least-squares estimator can be 
obtained by minimizing the sum of squared errors 


J N R(0)=T{z(j)-g[x(jXO]} z 

7=1 

(19) 

To demonstrate this estimation method the out-of-phase component, C N ^ , shown in Fig. 2, will be used as an 
example. In this case the regression equation is 



U) = c Nq 


(oo) - a 


i+A 2 (/) 


( 20 ) 


To check the regression, parameter estimates can be used to compute an estimated C N using Eq. (12). 

Measured C N and estimated C N are plotted in Fig. 3. These results show a good fit to the measured data. 

For a demonstration of model prediction capability and model validation, Fig. 4 shows the variation of C N 
with a at/= 0.86 Hz, a A - 10°, and two nominal values of a 0 = (14°, 26°). In both cases the shape of the 
predicted data forms a regular ellipse reflecting the underlying linear aerodynamic model structure. The 
measured C N shows some deviation from a regular ellipse, or linear behaviour, in the pre-stall and stall 
regions. 
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3.3 Two-Step Linear Regression 

An aerodynamic model of an aircraft performing a one degree-of-freedom oscillatory motion about one of its 
body axes can be formulated in terms of the in-phase and out-of-phase components, Ref. [12]. Using subscript 
a to represent the appropriate force or moment, the model can be written for rolling motion as 


for pitching motion as 

Cap = Cap (°°) sin a - af x sin a 
C.. =C a (oo)-a/ 0 sina 

p A ’ (21a) 

and for yawing motion as 

Ca a ~ C °a (°°) _ Cfl 

C aq =C (oo)-a/ 0 

9 q (21b) 

where 

Cap = c ap (co)cos a - afi cos a 
C n -C n (cc) + afocosa 

° r rK ’ (21c) 

With the expression, 

2, 2 
r, k 

h i+r y 

/o “ 1 + r,V < 22 > 

tfk 2 _ 1 

1 + T x k 2 1 + T x k 2 


(23) 

equations (21a), (21b) or (21c) can be rearranged into a set of equations for m different values of k as 


where for rolling oscillations 

y{j) = %+a l x(j), j = l,2,...,m 

(24) 
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x = C ap> y = C a p 

%= c a p (») + a 1 (a + C^(co))sina 

a x = -r x (25) 

and similarly for the pitching and yawing oscillations. In the first step a linear regression is used in estimation 
of parameters < 2 0 and a x \ n Eq. (25) from measured in-phase and out-of-phase components at m different 
values of k, m > 2. 

The second step of regression follows from equations (21a) and (22) replacing 2 ^ by its estimated value. The 
resulting equations are 


yi(j) = d o + d i x iU)’ j = h 2 ,...,m 

y2U) = c o +d i x 2(j)’ j = l2,...,m 


where for rolling oscillations these terms are 


(26) 


yi(j) = c afl (j), y 2 (j) = c ap 0 ) 

x i ( j ) = -/1 (;) sina ’ x 2 U) = ~foU) sina 

d 0= C a fi ( CO ) sina ’ d \ = a sin «> C 0 = C a p (°°) (27) 

Similar expressions are obtained for pitching and yawing oscillations. More discussion on the development of 
regression equations and estimator properties can be found in Ref. [12]. 

Two-Step Regression was applied to roll oscillatory data from NASA Langley 14x22 Wind Tunnel 
experiment using an 18% scale F-16XL aircraft. The results of harmonic analysis are plotted in Fig. 5. Both 
in-phase and out-of-phase components indicate frequency dependence for 25° < a 0 < 45°. The coefficient of 
determination confirms linearity of the aerodynamic model over the set of a 0 considered with the exception of 
a 0 = 45°. The linear dependence of both components is shown in Fig. 6. The slope of these data is equal to the 
time constant Tj. All parameter estimates for a 0 = 34° are presented in Table I. Variation of these parameters 
with a is given in Fig. 7, and a comparison of measured and predicted rolling moment coefficient is shown in 
Fig. 8. 


Table I. Two-Step Regression Parameter estimates for unsteady model at ao = 34° and oia - 10°. 


Step 1 

Step 2 

Tl 

a 

Cip 

Cip 

12.8 

0.331 

-0.099 

-0.039 

(0.22) 

(0.0041) 

(0.0074) 

(0.0084) 


3.4 Equation Error Method 

The Equation Error (EE) Method is in principle a linear regression. In general, it can be developed either in 
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the time domain or the frequency domain. However, considering the state-space representation in Eq. (4), the 
time domain approach is not possible because 77 and 77 are not measureable variables. 

For the frequency domain approach the model is given by Eq. (5) after expressing s as ico. In this case the 
model has the form 


C N (co) — 


-Aco 2 + C + Bico , . 

a(co) 

b x + ico 


(28a) 


z(j) = C N (j) + v( j\ j = 1, 2, ..., m 


(28b) 

where C N (co) and a(co) are the Fourier transforms of C N (t) and a(t ) , v(j) is the measurement noise, m is the 

number of frequencies at which the transformed input output data are known, and co is the angular frequency. 
In the Equation Error formulation unknown parameters are estimated from minimization of 


J EE (0) - ^ OX^l + 1(0 j ) + (Aa> 
7=1 


C -Bico ■ 


)«(7)| 2 , j = 1 , 2 ,..., m 


(29) 


In this formulation v(j) are residuals which encompass measurement noise and equation errors. Parameters 
A, B , and C are related to the aerodynamic coefficients by Eq. (6). 


An example of this method, Ref. [6], was applied to measurements obtained from pitch oscillation tests of a 
2.5% scale F-16XL model in a water tunnel. Also an alternative to the conventional single-frequency forced- 
oscillation test input was introduced in Ref. [6] by using Schroeder sweeps. This approach was found to 
provide substantially more efficient and more effective input for identification of unsteady models and was 
validated in Ref. [7]. Practical identification using a single-frequency input would require 6 runs, each at a 
different frequency, to ensure the four parameter model is adequately identified. For this example, a Schroeder 
sweep is used to allow model identification in one forced-oscillation run at a 0 - 42.5°. Figure 9 shows the 
input, a(t ) , and output, C N (t ) , displacements from their starting values at a 0 = 42.5°. The power spectrum 
of the Schroeder sweep is shown in Fig. 10. Schroeder sweeps provide a flat power spectrum over a specified 
frequency band with low peak-to-peak amplitudes. Amplitude of the Schroeder sweep is controlled by 
properly phasing each harmonic. Harmonic analysis provided the in-phase and out-of-phase components 
shown in Fig. 11. In this case, both components present strong frequency dependence for 30° < a < 70°, 
indicating an unsteady model is needed for that angle of attack region. Final estimates from the EE approach 
are shown with the Output Error (OE) results in the next section. In this study the EE estimates were used as 
initial estimates or starting values for the optimization process inherent in the OE method. 


3.5 Output Error Method (Frequency Domain) 

The output error method in the frequency domain can be applied directly using Eq. (28). In this case 
parameter estimates are obtained by minimization of the mean square output error. 
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m # 

Joe(0)= Zv (j)v(j) 
j = 1 

(30) 

where 


v(j) = z(j)- 


-Aco 2 ( j) + C + Bioj(j) 
b\ + ico(j) 


(31) 


The OE approach was applied to the same wide -band forced-oscillation data as the EE method discussed 
previously. Figure 12 shows the parameter estimates and their standard errors over the range of angle of attack 
where unsteady behaviours were observed. These estimates are consistent with the values obtained from the 
EE method. Estimate from the two methods in the frequency domain (f.d.) are compared in Table II for a 0 = 
42.5°. 


Table II. Parameter estimates for unsteady model at ao = 42.5°. 



Equation Error (f.d.) 

Output Error (f.d.) 

Output Error (t.d.) 

e 

6 

s{6) 

0 

s{6) 

0 

s{6) 

A 

0.939 

0.063 

0.895 

0.082 

not applicable 

B 

2.571 

0.055 

2.521 

0.064 

C 

-0.090 

0.051 

-0.101 

0.019 

bi 

0.138 

0.026 

0.144 

0.013 

Cncc 

-0.65 

0.39 

-0.70 

0.15 

- 0.35 

0.025 

CNq 

2.79 

0.19 

2.66 

0.24 

2.60 

0.045 

a 

-3.09 

0.39 

-3.09 

0.16 

-2.80 

0.025 

Ti 

21.53 

4.11 

20.57 

1.90 

17.48 

0.29 


Both the parameters directly estimated for the transfer function given by Eq. (28a) and those indirectly 
computed using Eq. (6) are shown in Table II. The two methods show very good agreement for the parameter 
estimates and the output error method obtained slightly better standard errors. Validation of these estimates is 
demonstrated by application to test data not used for estimation. Fig. 13 shows the output error estimates 
applied to ramp-and-hold data for a model validation test. The time histories show a very good match. 

3.6 Output Error Method (Time Domain) 

The estimation methods presented so far have been applied to linear aerodynamic models. These methods can 
be used to estimate all damping and cross derivative terms, however, an OE estimation method is needed 
when more general nonlinear aerodynamic model structures are required. In Ref. [14], a general mathematical 
model structure and modelling methodology was presented to address the nonlinear unsteady case. The 
methodology suggests using the least complex model that provides an adequate representation of the 
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aerodynamic response. With this approach, a progression is made toward more complex models and 
correspondingly more complex experiment design. The simplifying assumptions associated with the least 
complex model are removed in progression, only as required to achieve model adequacy. In general, this 
approach should lead to the best predictive model. 

An example from Ref. [14] using measurements from a single-axis forced-oscillation experiment in roll with a 
2.5% scale F-16XL model is presented here to demonstrate the method. An additional example using a 
transport configuration during single-axis yaw oscillations can be found in Ref. [13]. For the F-16XL 
example, measurements were obtained at six non-dimensional frequencies, k = [0.066, 0.095, 0.131, 0.160, 
0.197, 0.262], for angle of attack of 37.5°, during large amplitude (^ A =30°) roll oscillation tests. As shown in 
previous examples harmonic analysis provides key diagnostic information on the adequacy of a linear 
aerodynamic model. In this example some insight is also provided on the degree of nonlinear behaviour. In 
Fig. 14 the top plot shows the harmonic analysis for this case over a large range of angle of attack. Frequency 
dependent behaviour is observed for 25° < a < 40°, indicating an unsteady model is needed for that angle of 
attack region. The middle plot presents R 2 for a first harmonic model (the linear case). This plot indicates that 
for the lower frequencies the linear model does not completely explain the total variation present in the data. 
The lower plot presents R 2 for 1 st , 2 nd , and 3 rd order harmonic models. A dramatic improvement in R 2 occurs 
only when a 3 rd order harmonic model is used. For the k = 0.066 case, R 2 changes from 0.78, with the 1 st 
harmonic model, to 0.97 when a 3 rd order harmonic model is used. One contribution to the nonlinear 
behaviour is from the static terms in the model. Fig. 15 shows the fairly severe nonlinear character of the roll 
moment over angle of attack and sideslip. The surface is relatively smooth and linear at low and very high 
angles of attack; however in the mid angle-of-attack range, 20° < a < 50°, fairly severe static nonlinearities 
occur. 

As discussed in Ref. [14], for a single-degree-of-freedom roll test it is reasonable to assume that C x can be 
expressed as C x (/?, p) and that a slightly more general form than Eq. (4) can be expressed as 


f,(t) = -b 1 (j3)r 1 (t)-a(j3)j3 

C a (0 = C a (go; P) + l -C (oo; P)p{t) + 7/(0 

v ' (32) 


A broad class of nonlinear unsteady responses can be modelled with this basic model structure where the 
static term, steady-flow damping term, and indicial term parameters can vary nonlinearly with /? . In this 
example, it is assumed that each of the four unknown parameters (C z ( oo;/?), C Xp (? o;/?), a({3), /?(/?)) can be 
represented by polynomials in (3 . In general, experiments can be designed to identify models for the static, 
damping, and indicial terms separately, Ref. [14], although few facilities are capable of executing the motions 
required. For the OE formulation in the time domain, unknown parameters are estimated through 
minimization of 


j OE (V)=iv 2 (i) 

i=\ 


i 


where v(i) = C a (i)-C a (co-,/3(i))--C ap (co-,/3(i))p(i)-tj(i) (33) 

and r/(i) is computed from the state equation in Eq. (32). In this example only conventional single-frequency, 
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single-axis, forced-oscillation test data is used. Figure 16 shows sample time histories of the large amplitude 
oscillatory data at a 0 = 37.5° and J3 0 = 0°. The figure shows four cycles of oscillation for four frequencies, k = 
[0.066, 0.095, 0.131, 0.160]. Each frequency of the input and output measurements are stacked to ensure all 
frequency content is included in the analysis. Integration of Eq. (30) must be performed respecting the 
different conditions at each frequency. The first cycle of predicted response is eliminated to remove any initial 
transients and ensure that steady harmonic motion is used in calculating residuals. The amplitude, (f) A , is the 
same for each frequency; however, the resulting a and ft oscillations change during roll oscillations according 
to the kinematic relationships: 


P = sin 1 {sin(^)sin(# 0 )cos(^ 0 ) -cos(^)sin(^ 0 )} 

* = tan- 1 |cos W tan W + si "W si "»o> } 

[ cos(<9 0 )cos(^ 0 ) j ( 34 ) 


where 9 0 and y/ 0 are the nominal pitch and yaw offsets, respectively. The distortion of the rolling moment 
response from a steady harmonic wave, in particular for the low frequency cases, indicates nonlinear 
behaviour. 


Under the constraint of planar harmonic data, an iterative 2-stage process shown in Fig. 17 can be used to 
identify the complete model. Before starting the process, initial parameter values for a linear model are 
estimated using one of the methods described previously. Small amplitude motions tend to be well represented 
by linear models. Polynomials of the static values should be well modelled since these values can be directly 
measured in the wind tunnel; however for this example, the static terms were estimated from the dynamic 
measurements. 


To start the 2-stage process, first model structure determination of the unsteady term polynomials, a(P) and 
b(P % is performed using stepwise regression (SR) described in Ref. [5]. SR is an extension of linear regression 
that includes identification of statistically significant mathematical model structures. To setup the regression 
equation, a variable y(t) is formed by subtracting the static term and approximated steady-flow damping term 
from the force or moment measurement, C a (t), as 


y(t) = c a (t ) - C a (co;/?) - ^C ap (00 -P)p{t) 

Then the regression equation, showing the discrete measured values at time t(i), is given as 


( 35 ) 


y E ( i ) = -\b, ((3 e ( i))y E ( i ) + a(j3 E ( i))j3 E (/)] + s n (0, i = l,2,...,N 

( 36 ) 

where index E indicates the measured values and £ r] (i) is an equation error. Derivatives of measurement, y E , 
are formed using a locally smoothed numerical derivative, Ref. [5]. Parameter estimates are then updated with 
output error estimation. State and measurement equations for OE estimation are 
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m = -by (/?)/0 - a{p)P E ; n(t = 0) = 7(0) 

y E (t) = C aE (i) - C a mp E (i)) -Lr {co- fi E {i))p E (f) + v(i), i = 1,2 

V p (37) 

A second stage identification is done since C a ^ (oo; /?) is only known approximately and the estimates of a(/3) 

and b(P ) are based on this a priori information. To update the damping term, a second variable is formed 

z(t) = C a (t)-C a (^j3)-? 7 (t) 

=^c, mm P(t) 

v ' ( 38 ) 

If the measured values and new estimates for the unsteady term are substituted into Eq. (38), a new regression 
equation is obtained as 


z E (i) = y C a p C 00 ;/?/?')/ E (i) + S z (i), i = l,2,...,N 


( 39 ) 


The model structure and parameters for 

the steady-flow damping term can be Tab | e m . Nonlinear model for C, at « 0 = 37.5°, /3 0 = 0°. 

estimated using Eq. (39). This 
completes the first iteration of the two 
stage identification. In the next 
iteration, a(fi) and b(J3) can be 
estimated again using the new values 
and model structure for the steady-flow 
damping term. This process is 
continued until parameter estimates 
have converged. 

Since harmonic analysis inferred that at 
least a cubic nonlinearity is present and 
the static component is a relatively 
large component in this case, a cubic 
representation of the static data was 
initially implemented. It was found, however, that a 5 th order polynomial was required for an adequate static 
model. Table III shows the final estimated model from the 2-stage process applied at a 0 = 37.5° and the 
corresponding overall R 2 . The parameter covariance matrix for this case contained some large pair-wise 
correlations indicating some identifiability issues. Figure 18 shows the measured and predicted rolling 
moments for the four frequencies considered. These graphs all show distortion from the regular elliptic shape 
associated with linear response. For each frequency, the nonlinear model prediction matches the response 
well. At the lowest frequency, k = 0.066, the distortion or nonlinearity is greatest, but the model match is still 
very good. Although the fit is very good the model in Table III is not completely satisfactory due to some 
inconsistency with small-amplitude linear analysis. Ideally the nonlinear model should represent the linear 
case as a subset then, for example, at (3 - 0, the nonlinear model estimate, C lp (oo; (3 = 0) , should be close to the 


Parameter 

Nonlinear Model 

C,(oo; P) 

0.007 +0.317/? +0.202^-12.786^ 
-2.371^+91.123/ 

<V°°;j0) 

-0.786 -1.791/ +97.616/ 
(0.0045) (0.34) (4.48) 

a(P) 

0.527 -43.668/ -507.228/ 
(0.0051) (0.39) (5.096) 

W 

0.624 -8.109/ 
(0.011) (0.14) 

II 

o 

3.63 

(0.027) 

R z 

0.97 
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corresponding linear model analysis of the small amplitude data, a value of -0.13. The differences may reflect 
a number of possible sources of error besides the model structure for the damping and unsteady terms, such as 
data accuracy and information content of the data. 


4.0 CONCLUDING REMARKS 

This paper presented several methods for estimating mathematical models useful for stability and control 
analysis that can be applied to both CFD simulations and wind tunnel measurements. These methods can 
accommodate general aerodynamic behaviors and characterize model parameter uncertainty. Uncertainty 
information significantly improves model validation tests by providing more information than conventional 
time history comparisons. Conventional stability and control derivatives are retained in the recommended 
model structures, even for the nonlinear case, to take advantage of the experience and significant knowledge 
base built up through years of aerospace engineering practice. 

Collaboration between computational and experimental researchers offers a significant opportunity to advance 
the technology in both disciplines. A number of modelling research areas present opportunity for collaboration 
that can be further developed. Dynamic test techniques in wind tunnels are limited to current facility 
capabilities and most facilities only allow planar harmonic motion. Using CFD simulations offers an 
opportunity to further the design of both dynamic inputs for experiments and desirable facility capabilities. 
Efficacy of simple improvements such as using wide-band inputs, oscillatory coning, and plunging dynamics 
can be relatively quickly and inexpensively assessed with CFD experiments. Some flight conditions are not 
easily reproduced in a tunnel and some conditions are either too transient or too dangerous for flight test such 
as with the study of transport loss of control problems. Collaboration between CFD and SID can facilitate 
both discovery and analysis of these difficult-to-model flight conditions. This approach can lower concerns for 
pilot and vehicle safety and potentially reduce costs. 
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S = 5.90 ft 2 , c = 0.915 ft, b = 6.85 ft 

Figure la. Model geometry for NASA Figure lb. Three-view drawing of 2.5% F-16XL water- 

GTM, experimental sub-scale aircraft. tunnel model. 





Figure 2. Harmonic analysis for normal-force coefficient, Figure 3. Results of nonlinear regression for measured 

NASA GTM configuration, a A = 10°. out-of-phase components of NASA GTM, ag= 18°, 

04=10°. 
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Figure 4. Comparison of measured and predicted 
normal-force coefficient, NASA GTM, f = 0.86 Hz, 
c^=10°. 
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Figure 5. Harmonic analysis for roll-moment coefficient, 
18% F-16XL configuration, a A = 10°. 



Figure 6. Measured vs Estimated roll- 
moment coefficients, 18% F-16XL, a A =10°. 
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Figure 7a. Two-Step Regression estimates, 18% F-16XL, 
0 ^= 10 °. 



Figure 7b. Two-Step Regression estimates, 18% F-16XL, 
a A =10°. 



Figure 7c. Two-Step Regression estimates, 18% F-16XL, 
0 ^ 4 = 10 °. 
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Figure 8. Two-Step Regression model of roll-moment predicted and measured data, 18% F-16XL, 0^=10° 



Figure 9. Measurements of a perturbations from 
Oo = 42.5° and normal force during wide-band 
experiment with 2.5% F-16XL in water tunnel. 


Figure 10. Harmonic content of a wide-band 
input, 2.5% F-16XL in water tunnel, a A = 5°. 
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Figure 12a. Normal Force model parameters and 
2-a bounds, using OE, 2.5% F-16XL, a A = 5°. 
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Figure 13. Model validation of C N response to 
ramp (cfy=40°-45°), non-dim pitch rate of 0.03, 
2.5% F-16XL. 
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Figure 15. Surface plot of roll moment static 
measurements, 2.5% F-16XL. 
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Figure 14. Harmonic Analysis of 2.5% F-16XL 
during large amplitude roll oscillations, (j) A = 30°. 
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Figure 16. Time histories of the large amplitude oscillatory 
data at cco = 37.5°, p 0 - 0°, and four frequencies, k = [0.066, 
0.095, 0.131, 0.1601, 2.5% F-16XL. 
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Figure 17. Block diagram of model identification procedure using 
stepwise regression (SR) and an output error (OE) estimation. 
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Figure 18. Comparison of measured and computed rolling- 
moment coefficient from time domain OE method at ao = 37.5°. 
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